\chapter{Descrizione del codice}

\includegraphics[width=0.3\columnwidth]{logo}

La libreria in linguaggio C++ \texttt{Heat-Spectral} è stata progettata per risolvere problemi parabolici bidimensionali, su griglie strutturate cartesiane, attraverso il metodo spettrale SEM-NI.
Il progetto è stato sviluppato a partire da un codice già esistente, \texttt{MPI-Spectral}, che implementa il metodo SEM-NI, ma per problemi ellittici e con condizioni al bordo solo di tipo Dirichlet, secondo un'ottica multidominio.
I principali sviluppi di cui ci siamo occupati hanno riguardato:
\begin{itemize}
	\item l'estrazione del solutore ellittico monodominio, indispensabile per poi procedere con ulteriori ampliamenti;
	\item l'implementazione delle condizioni di tipo Neumann non omogenee, in modo da riuscire ad affrontare anche simulazioni di problemi più complessi;
	\item l'estensione al caso di problemi parabolici, grazie all'implementazione del $\vartheta$-metodo;
	\item l'implementazione di un'interfaccia grafica, attraverso la generazione di opportuni script, che consente di generare immagini e video rappresentanti l'evolversi della soluzione numerica.
\end{itemize}
Il progetto è stato messo fin dall'inizio sotto controllo di versione grazie al sistema \texttt{subversion}.
Questa pratica consente di registrare su server ogni modifica del codice in fase di sviluppo, in modo da poter ripristinare vecchie versioni in caso di errori, oltre a permettere una migliore interazione tra i vari sviluppatori, che possono lavorare contemporaneamente allo stesso progetto senza rischiare di interferire tra loro.
Il server che è stato scelto per ospitare il nostro progetto è la piattaforma di Google Code, che fornisce, a titolo gratuito, lo spazio e le funzionalità per sviluppare strumenti software \emph{open source}.
Grazie a questo supporto, abbiamo avuto a disposizione una pagina web liberamente accessibile, dalla quale chiunque può ottenere tutte le informazioni relative al nostro progetto, oltre ad avere la possibilità di scaricarne i sorgenti e la relativa documentazione\footnote{\url{http://code.google.com/p/heat-spectral/}}.\\
A tal proposito è bene sottolineare che, per produrre un manuale rigoroso e coerente, ci siamo avvalsi del generatore automatico di documentazione \texttt{doxygen}.
Questa applicazione consente di produrre una ricca guida all'utilizzo della libreria, semplicemente a partire dai commenti inseriti all'interno del codice sorgente.
Più avanti discuteremo nel dettaglio gli aspetti tecnici degli strumenti software esterni utilizzati.\\
In effetti, oltre a quelli appena citati, ci siamo avvalsi di altri supporti, più strettamente legati al funzionamento della libreria stessa.
Per esempio, per quanto riguarda la risoluzione del sistema lineare associato al problema differenziale, abbiamo utilizzato la libreria \texttt{UMFPACK}~\cite{UMFPACK}, che effettua una fattorizzazione LU in formato sparso \texttt{CSC} (\emph{Compressed Sparse Column}), rendendo efficiente l'inversione delle matrici.
Per ulteriori dettagli sulle tecniche di memorizzazione di matrici sparse si può approfondire consultando \cite{FORMAGGIA}.\\
L'interfaccia grafica, invece, è stata implementata attraverso la generazione di script che permettano chiamate immediate al programma di visualizzazione grafica \texttt{gnuplot}~\cite{GNUPLOT}.
Di fatto sono stati costruiti dei metodi che, grazie alla possibilità di conoscere tutte le informazioni relative al problema in fase di risoluzione, sono in grado di costruire script personalizzati i quali, con una sola chiamata, permettono di generare immagini della soluzione a tempi fissati.
In realtà, sfruttando le potenzialità del software per la generazioni di filmati \texttt{MEncoder}~\cite{MPLAYER}, è anche possibile ottenere, in maniera piuttosto rapida, un video relativo all'intero evolversi del problema parabolico.
In effetti, data la caratteristica evolutiva dei sistemi che intendiamo simulare, è naturale voler visualizzare l'andamento temporale con una sequenza che colga ogni singolo istante.\\
Infine, per quanto riguarda l'aspetto maggiormente legato alla tecnica di programmazione, si è cercato di templetizzare al massimo la struttura del codice, dove possibile, in modo da garantire sufficiente generalità e astrazione.
A tal proposito, sono state di enorme aiuto le accurate descrizioni della sintassi del C++ contenute in~\cite{YANG}.
Inoltre, al fine di rendere più pulito l'utilizzo della libreria, si è lasciata all'utente la possibilità di intervenire solo sulle parti esterne del codice, dedicate alla costruzione del problema, rendendo privati i solutori principali.
In realtà, già il codice da cui è cominciato il nostro sviluppo era fortemente templetizzato e, pertanto, è stato necessario attenersi allo stile di programmazione preposto, in modo da garantire una certa omogeneità nel \emph{modus operandi}.
In particolare, il fatto di espandere una libreria già implementata, invece di cominciare il lavoro dal principio, ha comportato una difficoltà maggiore, poiché si è presentata la necessità di comprendere in ogni dettaglio il codice pre-esistente che, in questo caso, si è rivelato essere decisamente sofisticato.\\
Dopo questa breve panoramica, passiamo ora a mostrare le principali classi del codice, cercando di seguire un ordine di complessità crescente.


\section{Le classi di base}

Il nucleo centrale del codice è costituito dalla classe \texttt{Problem}, la quale contiene tutti i metodi per la costruzione delle matrici di massa e di stiffness, oltre alle funzioni necessarie alla loro inversione, che avviene sfruttando gli strumenti forniti dalla libreria \texttt{UMFPACK}.
È proprio la classe \texttt{Problem} a chiamare tutte le altre che compongono la libreria, in modo da impostare e risolvere il problema evolutivo che intendiamo sottoporle.
Vediamo allora, più nel dettaglio, quali sono le classi che vengono chiamate a svolgere tutte le operazioni necessarie.

\subsection{La classe \texttt{Point}}

Tra le classi geometriche, la classe \texttt{Point} è la più semplice: di fatto contiene soltanto la definizione di punto in un piano cartesiano.
Il tipo per le coordinate è templetizzato, in modo che sia possibile definire dei punti con delle coordinate intere (nel caso di griglia discreta), oppure con delle coordinate reali (nel caso di griglia continua).
Oltre a questo, sono naturalmente presenti dei metodi classici, in grado di svolgere tutte le operazioni relative alla gestione del punto, all'interno della dinamica complessiva della libreria.

\subsection{La classe \texttt{Poly}}

La classe \texttt{Poly} permette di rappresentare i polinomi, memorizzandone il grado e i coefficienti. Implementa anche le principali funzioni che coinvolgono i polinomi, come il calcolo degli zeri attraverso l'algoritmo di Newton-Horner (per i cui dettagli rimandiamo a \cite{NUMERICA}) e la generazione dei polinomi di Legendre necessari alle costruzione delle formule di integrazione numerica, come descritto nel precedente capitolo.
Infine, implementa anche le operazioni di integrazione e differenziazione di polinomi, procedimenti indispensabili alla risoluzione del problema differenziale attraverso il metodo spettrale.

\subsection{La classe \texttt{TensorProd}}

La classe \texttt{TensorProd} implementa il prodotto tensoriale tra due polinomi monodimensionali, del tipo:
\[
	\left(f \otimes g\right)(x,y) = f(x) g(y).
\]
In questo modo è possibile costruire la base spettrale bidimensionale, definita come prodotto tensore tra polinomi costruiti sulla griglia GLL.
È naturale che questa classe non viene utilizzata direttamente dall'utente che si interfaccia all'utilizzo della libreria, ma è di fondamentale importanza nella gestione delle integrazioni numeriche sulla griglia bidimensionale.

\subsection{La classe \texttt{Integrator}}

La classe \texttt{Integrator} implementa le formule di integrazione gaussiana sui nodi GLL della griglia rettangolare di riferimento.
Tra i principali attributi presi in ingresso ci sono il numero di nodi utilizzati dalla formula di quadratura in ogni elemento e il numero massimo di iterazioni del metodo di Newton per il calcolo delle radici del polinomio di Legendre.\\
Questa classe non può essere utilizzata direttamente dall'utente, dato che sarebbe necessario conoscere la struttura della mesh.
Inoltre, per avere l'integrazione su uno specifico elemento, è necessario applicare la trasformazione che mappa l'elemento di riferimento in quello specifico, come in (\ref{eq:gll}).
Questi due problemi sono risolti, rispettivamente, dalle classi \texttt{Mesh} e \texttt{Map2d}.

\subsection{La classe \texttt{Map2d}}

La classe \texttt{Map2d} implementa la trasformazione $\boldsymbol{\varphi}$ tra l'elemento di riferimento $\hat{\varOmega} = [-1,1]\times[-1,1]$ e un qualsiasi elemento rettangolare assegnato, il quale è caratterizzato da $2$ dei suoi vertici: quello in basso a sinistra e quello in alto a destra.
In particolare, la trasformazione diretta mappa l'elemento di riferimento in quello assegnato mentre, per quanto riguarda l'inversa, vale il contrario.
Infine sottolineiamo come, anche in questo caso, la tipologia di mappa che intendiamo applicare, se diretta o inversa, è passata al costruttore per mezzo di un template, sempre secondo la filosofia di programmazione alla quale accennavamo in precedenza.

\subsection{Le classi \texttt{Functor} e \texttt{Product}}

La classe \texttt{Functor} permette di effettuare la composizione tra una funzione $f$ e una mappa affine $\boldsymbol{\varphi}$:
\[
	f\left(\varphi_1(x,y),\varphi_2(x,y)\right),
\]
in modo da poter valutare la funzione su uno specifico elemento, anziché esclusivamente su quello di riferimento.
Naturalmente non c'è alcun bisogno di sottolineare la fondamentale importanza di questa classe nell'ottica globale del codice.\\
Invece, per quanto riguarda la classe \texttt{Product}, questa permette di effettuare il prodotto tra due funzioni $f(x,y)$ e $g(x,y)$, indispensabile poi per costruire il prodotto scalare nello spazio $L^2(\varOmega)$.


\section{La classe \texttt{Mesh}}

La classe \texttt{Mesh} permette di memorizzare la griglia di calcolo bidimensionale su un dominio rettangolare per il metodo SEM-NI, e fornisce i principali strumenti di calcolo utilizzabili nella fase della costruzione del problema da parte dell'utente.\\
Tra gli attributi privati, mostriamo quelli più importanti:
\begin{description}[\texttt{Point$<$double$>$ \_A}:]
	\item[\texttt{int \_nx}:] numero di elementi sull'asse x;
	\item[\texttt{int \_ny}:] numero di elementi sull'asse y;
	\item[\texttt{int \_nnodes}:] numero di nodi utilizzati dalla formula di quadratura GLL;
	\item[\texttt{Point$<$double$>$ \_A}:] punto del dominio in basso a sinistra; 
	\item[\texttt{Point$<$double$>$ \_B}:] punto del dominio in alto a destra;
	\item[\texttt{double \_legDelta}:] minimo incremento nella formula di Newton utilizzata per ottenere i nodi GLL;
	\item[\texttt{double \_legEps}:] residuo critico nella formula di Newton utilizzata per ottenere i nodi GLL;
	\item[\texttt{double *\_vertexX}:] vertice sull'asse x;
	\item[\texttt{double *\_vertexY}:] vertice sull'asse y.
\end{description}
Sono presenti anche numerosi metodi pubblici, la maggior parte dei quali dedicata alla gestione delle principali proprietà geometriche:
\begin{itemize}
	\item ritornare la numerazione dei nodi, sia a livello locale che a livello globale, e la numerazione degli elementi;
	\item ritornare le coordinate dei nodi, sempre sia a livello locale che a livello globale;
	\item ritornare la restrizione di una certa funzione o di un vettore su un elemento assegnato;
	\item ritornare la traccia al bordo di una funzione;
	\item ritornare la base lagrangiana;
	\item interpolare una funzione $f$ sulla mesh;
	\item ricavare le derivate spaziali di oggetti di tipo \texttt{TensorProd} oppure \texttt{Solution};
	\item calcolare il prodotto scalare $L^2$ tra funzioni.
\end{itemize}
Infine, è anche presente un metodo, chiamato \texttt{save\_to\_file}, che si occupa  di salvare in un file di testo i valori di un vettore, in corrispondenza dei rispettivi nodi della mesh.
Di fatto questo metodo permette di salvare il vettore soluzione in un formato utile per una successiva visualizzazione grafica.


\section{La classe \texttt{UMFMatrix}}

La classe \texttt{UMFMatrix} permette di memorizzare le matrici in formato sparso \texttt{CSC}, costruendo, al posto di un'intera matrice, soltanto tre vettori:
\begin{itemize}
	\item il vettore \texttt{\_Ax} dei valori non-nulli nella matrice;
	\item il vettore \texttt{\_Ap} degli indici del primo elemento di ogni colonna;
	\item il vettore \texttt{\_Ai} degli indici della posizione verticale di ogni elemento.
\end{itemize}
Lo scopo principale di questa classe è di interagire con la libreria \texttt{UMFPACK} allo scopo di fattorizzare le matrici e risolvere sistemi lineari del tipo:
\begin{displaymath}
	\A\vec{x} = \vec{b}
\end{displaymath}
Passiamo dunque a mostrare i metodi principali che si occupano di svolgere queste operazioni.
\begin{description}
	\item[\texttt{multiply}:] Riportiamo il segmento di codice in modo da mostrarne il funzionamento.
		\begin{c++}
void UMFMatrix::multiply(double const * V, double * AV) const {
    for(int i=0;i<_sizeA;++i) {
        AV[i]=0.0;
    }
    for(int i=0;i<_sizeA;++i) {
        for(int j=_Ap[i];j<_Ap[i+1];++j) {
            assert(j<_sizeNZ);
            assert(_Ai[j]<_sizeA);

            AV[i]+=_Ax[j]*V[_Ai[j]];
        }
    }
}
		\end{c++}
	Il metodo permette di moltiplicare la matrice $\A$ per il vettore $\mathrm{V}$, memorizzando il risultato nel vettore $\A\mathrm{V}$;
	\item[\texttt{factorize}:] Questo metodo non prende in ingresso alcun parametro:
		\begin{c++}
void UMFMatrix::factorize() 
		\end{c++}
		ed è in grado di fattorizzare la matrice sfruttando i metodi di \texttt{UMFPACK};
	\item[\texttt{solve}:] Questo metodo prende in ingresso il vettore $b$ dei termini noti e il vettore $x$ nel quale, poi, andrà a scrivere la soluzione:
		\begin{c++}
void UMFMatrix::solve(double const * b, double * x) const
		\end{c++}
		La risoluzione del sistema lineare, secondo quanto già affermato, è effettuata attraverso gli algoritmi ottimizzati della libreria \texttt{UMFPACK}.
\end{description}
Relativamente ai metodi secondari di questa classe, possiamo accennare a quelli che restituiscono gli elementi dei tre vettori in cui viene memorizzata la matrice, oltre a quelli che permettono di ridimensionare dinamicamente le lunghezze di tali vettori.\\
Per ultimo, allo scopo di ottenere uno stile di programmazione più pulito, è anche stato effettutato l'\emph{overloading} dei principali operatori relativi alle strutture algebriche.
A questo proposito ricordiamo quello relativo alle operazioni di copia e di somma tra matrici, oltre all'indispensabile operazione di prodotto di una matrice per uno scalare.\\
Per finire, abbiamo anche pensato di introdurre una semplice funzione in grado di riportare in un file di testo il \emph{pattern} della matrice.
In questo modo, grazie a un qualsiasi programma di visualizzazione grafica, oppure importando il file in un software di calcolo scientifico, è immediato rappresentare la struttura sparsa della matrice stessa.
È grazie a questo metodo che siamo in grado di evidenziare le differenze nella struttura della matrice al variare del numero di elementi spettrali.


\section{La classe \texttt{Problem}}

La classe \texttt{Problem} è, come anticipato in precedenza, quella più importante di tutta la nostra libreria, poiché permette all'utente di definire e risolvere il problema evolutivo.
Cominciamo col mostrare la firma del costruttore:
\begin{c++}
template<class BIL,class F>
class Problem {
    public:
        Problem(Point<double> const & AP,
                Point<double> const & BP,
                int nx,
                int ny,
                int nnodes,
                double deltat,
                double theta,
                int nmax,
                double legDelta,
                double legEps,
                int order=1);
\end{c++}
Evidenziamo come le classi \texttt{BIL} e \texttt{F}, che portano l'informazione relativa a forma bilineare e funzionale, vengano passate a template, secondo lo schema ampiamente descritto in precedenza.
Il costruttore si occupa di inizializzare i parametri necessari alla costruzione della mesh e, quindi, delle matrici.\\
Ma veniamo a mostrare i metodi principali di questa classe.
\begin{itemize}
	\item Il metodo \texttt{initialize}, che genera il problema il problema evolutivo e costruisce le strutture algebriche necessarie alla sua risoluzione.
	\item Il metodo \texttt{solve}, che risolve per ogni istante di tempo il problema parabolico.
	\item I metodi di \emph{post-processing}, che permettono il salvataggio in un file della soluzione a ogni passo temporale e la generazione degli script per la sua visualizzazione grafica.
\end{itemize}
Chiaramente, a loro volta, ciascuno di questi metodi principali richiama diversi metodi secondari, i quali sono i veri attori del processo risolutivo.
Prima di analizzarli nel dettaglio, però, torniamo per un attimo alla formulazione algebrica del metodo SEM-NI per problemi parabolici, quella descritta in (\ref{eq:thetametodo}).
Abbiamo infatti la necessità di mostrare come l'approccio analitico venga interpretato in chiave algoritmica.
Per semplicità, ipotizziamo inizialmente che il funzionale sia costante nel tempo:
\begin{equation}\label{eq:ipotesicost}
	\vec{f}^k = \vec{f}^{k+1}, \qquad \forall k>0.
\end{equation}
In questo modo, la (\ref{eq:thetametodo}) diventa:
\begin{equation}\label{eq:matricek}
	\underbrace{\left(\dfrac{\M}{\Delta t} + \vartheta \A \right)}_\K \vec{u}^{k+1} = \underbrace{\left(\dfrac{\M}{\Delta t} - (1-\vartheta)\A \right)}_\mathrm{G} \vec{u}^k + \vec{f}^k.
\end{equation}
Ora, se osserviamo che
\[
	\K^{-1}\mathrm{G} = \left[\dfrac{\M}{\Delta t} +\vartheta\A\right]^{-1} \left[\left(\dfrac{\M}{\Delta t}+\vartheta\A\right) - \A \right] = \K^{-1}(\K-\A) = \mathrm{I}-\K^{-1}\A,
\]
allora possiamo ottenere
\begin{eqnarray*}
	\vec{u}^{k+1}	&=	& (\mathrm{I}-\K^{-1}\A)\vec{u}^k + \K^{-1}\vec{f}^k\\[2mm]
				&=	& \vec{u}^k + \K^{-1}(\vec{f}^k-\A\vec{u}^k)\\[2mm]
				&=	& \vec{u}^k + \K^{-1}\vec{r}^k\\[2mm]
				&=	& \vec{u}^k + \boldsymbol{\delta}^k,
\end{eqnarray*}
dove $\vec{r}^k$ indica il vettore dei residui al passo $k$, cioè
\[
	\vec{r}^k = \vec{f}^k - \A\vec{u}^k,
\]
mentre $\boldsymbol{\delta}^k$ rappresenta l'incremento necessario per aggiornare la soluzione dal passo $k$ al passo $k+1$.\\
Il sistema lineare che dovrà essere risolto a ogni passo, pertanto, risulta essere
\begin{equation}\label{eq:sistema}
	\K \boldsymbol{\delta}^k = \mathbf{r}^k
\end{equation}
ma, con evidente soddisfazione, notiamo che la matrice $\K$ non deve essere ricalcolata a ogni istante temporale e, di conseguenza, può essere fattorizzata una volta per tutte rendendo la risoluzione del problema molto più immediata.
Incidentalmente osserviamo che l'incremento della soluzione si traduce, quindi, in una semplice somma vettoriale.\\
In realtà la nostra libreria è anche in grado di trattare il caso di forzanti tempo-dipendenti, in cui non è più valida l'ipotesi (\ref{eq:ipotesicost}).
In questo caso la formulazione algebrica diventa:
\begin{equation}
	\vec{u}^{k+1} = \vec{u}^k + \K^{-1}\left(\vec{r}^k + \vartheta\left(\vec{f}^{k+1}-\vec{f}^k\right)\right)
\end{equation}
e, di conseguenza, notiamo come ora il vettore incremento sia la soluzione del sistema lineare
\[
	\K\boldsymbol{\delta}^k = \vec{r}^k + \vartheta(\vec{f}^{k+1}-\vec{f}^k).
\]
Questa modifica ha comportato qualche complicazione a livello implementativo poiché, in questo caso, è necessario ricalcolare il vettore $\vec{f}$ ad ogni istante di tempo.
In realtà, però, tale vettore non contiene esclusivamente le informazioni relative alla forzante evolutiva, bensì anche quelle collegate all'imposizione delle condizioni al contorno che, al contrario, rimangono inalterate nel tempo.
In conclusione, la modifica a ogni passo del vettore $\vec{f}$ dovrà tener conto di entrambi i contributi e, a questo scopo, sono giustificati gli sforzi in fase di programmazione per ottenere un codice snello, che consenta di affrontare anche questo scenario. 

\subsection{Il metodo \texttt{initialize}}

Il metodo \texttt{initialize} si occupa di generare le strutture algebriche necessarie alla risoluzione del problema SEM-NI in versione algebrica.
Cominciamo con il mostrarne la firma.
\begin{c++}
void initialize(UMFMatrix & A,
                UMFMatrix & M,
                UMFMatrix & K,
                std::valarray<double> & b,
                std::valarray<double> & bold,
                std::valarray<double> & u,
                std::valarray<double> & uold,
                std::valarray<double> & residual,
                std::valarray<double> & delta,
                std::valarray<double> const & bc);
\end{c++}
Possiamo notare come figurino tutte le matrici e tutti i vettori che abbiamo introdotto nella precedente sezione.
In aggiunta, compare anche un vettore costante, che serve a immagazzinare l'informazione relativa alle condizioni al contorno, sia che si tratti di condizioni essenziali, sia quando abbiamo la necessità di imporre delle condizioni naturali.\\
Scendendo più nel dettaglio del funzionamento di \texttt{initialize}, mostriamo come nel corpo di questo metodo figurino le seguenti istruzioni:
\begin{c++}
	A.resize_sizeA(_sizeA);
	M.resize_sizeA(_sizeA);
	K.resize_sizeA(_sizeA);
	b.resize(_sizeA);
	bold.resize(_sizeA);
	u.resize(_sizeA);
	uold.resize(_sizeA);
	residual.resize(_sizeA);
	delta.resize(_sizeA);
	matrixGen(A);
	massMatrixGen(M);
	forceGen(b,0);
	patchMatrix(A);
	bcGen(bc,b);
	initGen(u);
	kappaMatrixGen(A,M,K);
\end{c++}
Le prime nove righe di comando si occupano di dimensionare correttamente le matrici e i vettori implicati:
\begin{description}[\texttt{residual}:]
	\item[\texttt{A}:] Matrice di stiffness;
	\item[\texttt{M}:] Matrice di massa;
	\item[\texttt{K}:] Matrice del $\vartheta$-metodo, definita nella (\ref{eq:matricek});
	\item[\texttt{b}:] Vettore dei termini noti al passo $k+1$;
	\item[\texttt{bold}:] Vettore dei termini noti al passo $k$;
	\item[\texttt{u}:] Vettore della soluzione al passo $k+1$;
	\item[\texttt{uold}:] Vettore della soluzione al passo $k$;
	\item[\texttt{residual}:] Vettore dei residui;
	\item[\texttt{delta}:] Vettore degli incrementi.
\end{description}
Invece, la definizione vera e propria del problema viene eseguita dai metodi successivi.
Analizziamoli più nel dettaglio.
\begin{description}[\texttt{kappaMatrixGen}:]
	\item[\texttt{matrixGen}:] Costruisce la matrice di stiffness, senza considerare le condizioni al contorno.
	Per ottenere la miglior efficienza computazionale possibile, in accordo con la struttura \texttt{CSC} della matrice, il metodo costruisce innanzitutto il pattern della matrice, e solo in un secondo momento assegna i valori numerici, valutando la forma bilineare sui nodi.
	\item[\texttt{massMatrixGen}:] È il metodo analogo a \texttt{matrixGen} che, però, si occupa di costruire la matrice di massa.
	Come nell'altro caso, la costruzione della matrice è suddivisa nelle due fasi di definizione del pattern e di valutazione dei valori numerici da inserire negli elementi non nulli della struttura algebrica.
	\item[\texttt{forceGen}:] Genera il vettore dei termini noti, valutando il funzionale della forzante sui nodi della griglia, sia nel caso di forzante costante sia in quello più generale di forzante evolutiva.
	Non vengono invece presi in considerazione i termini dovuti alle condizioni al contorno.
	Questa operazione è delegata al metodo \texttt{solve}.
	\item[\texttt{patchMatrix}:] Si occupa di rendere la matrice di stiffness consistente con le condizioni al contorno di tipo Dirichlet.
	In corrispondenza di questi nodi, viene assegnato un valore unitario all'elemento sulla diagonale principale, e valore nullo a tutti i restanti elementi di quella specifica riga.
	In questo modo, di fatto, stiamo imponendo in maniera forte che la soluzione in quel preciso punto assuma il valore imposto dalla condizione di Dirichlet.
	\item[\texttt{bcGen}:] Analogamente al metodo \texttt{patchMatrix}, si occupa di rendere consistente il vettore dei termini noti con le condizioni al contorno di tipo Dirichlet e Neumann.
	Per svolgere questo compito il metodo si appoggia alla funzione \texttt{edgeBc}, che lavora singolarmente sui quattro bordi.
	Nel caso di condizione Dirichlet, viene sostituito il valore fisico del bordo mentre, nel caso di condizione Neumann, viene aggiunto il prodotto scalare $L^2(\varGamma_N)$ tra la funzione che la definisce e la base lagrangiana di riferimento.
	\item[\texttt{initGen}:] Costruisce la soluzione al passo iniziale $k=0$, effettuando semplicemente un'interpolazione della funzione dei valori iniziali sui nodi della griglia.
	\item[\texttt{kappaMatrixGen}:] Costruisce la matrice $\K$ a partire dalle matrici $\A$ e $\M$, secondo l'operazione definita in (\ref{eq:matricek}), e poi ne effettua la fattorizzazione LU, secondo le modalità definite in precedenza.
\end{description}
Secondo l'approccio algebrico alla risoluzione, nonostante queste operazioni di inizializzazione abbiano un elevato onere computazionale, una volta effettuate la risoluzione di ogni passo temporale è resa decisamente più immediata, considerata la stazionarietà della matrice $\K$.

\subsection{Il metodo \texttt{solve}}

Il metodo \texttt{solve} si occupa esclusivamente della risoluzione del sistema lineare (\ref{eq:sistema}) a ogni passo temporale $k$.
Grazie all'aver svolto tutte le operazioni preliminari in un momento precedente, esso riesce a essere decisamente più snello del previsto, come si può osservare dal numero limitato di istruzioni.
È proprio questo metodo che, in fase di risoluzione di un ben preciso problema evolutivo, può essere richiamato dall'utente all'interno del suo programma.
Mostriamone allora il contenuto, illustrandone il suo elementare funzionamento.
\begin{c++}
template<class BIL,class F>
void Problem<BIL,F>::solve() {
    _uold = _u;
    _A.multiply(_uold, _residual);
    _residual = _b - _residual;
    _K.solve(_residual, _delta);
    _u = _uold + _delta;
}
\end{c++}
Per cominciare, la soluzione al passo precedente viene memorizzata nel vettore \texttt{\_uold}.
Dopodiché si calcola il prodotto fra la matrice di stiffness e la soluzione, in modo da ottenere il vettore dei residui.
A questo punto, il metodo \texttt{solve} della classe \texttt{UMFMatrix} si occupa di calcolare la soluzione \texttt{\_delta} del sistema lineare (\ref{eq:sistema}), una volta che gli abbiamo fornito la matrice e il termine noto.
Infine è sufficiente effettuare una semplice somma di vettori, allo scopo di ottenere la soluzione \texttt{\_u} al nuovo passo temporale.\\
Per quanto riguarda il caso di forzante non costante, invece, è necessario che il metodo \texttt{solve} prenda in ingresso anche l'istante temporale:
\begin{c++}
template<class BIL,class F>
void Problem<BIL,F>::solve(double const & t) {
    _uold = _u;
    _bold = _b;
    forceGen(_b, t);
    bcGen(_bc, _b);
    _A.multiply(_uold, _residual);
    _residual = _theta*_b + (1-_theta)*_bold - _residual;
    _K.solve(_residual, _delta);
    _u = _uold + _delta;
}
\end{c++}
In questo modo stiamo sfruttando le potenzialità del linguaggio di programmazione, che ci consente di generare più metodi con lo stesso nome, in modo che siano differenziati dalla loro firma.
Questo stratagemma ci pare decisamente funzionale per la costruzione di un codice pulito dal punto di vista stilistico ma, allo stesso tempo, capace di operare in un regime di una certa complessità.\\
Diversamente dal caso precedente, ora ad ogni istante temporale vengono valutati i termini forzanti e i termini di bordo dovuti alle condizioni al contorno.
Di conseguenza, come è naturale aspettarsi, la risoluzione di un problema con forzante evolutiva non può che richiedere uno sforzo computazionale superiore.

\subsection{I metodi di post-processing}

I metodi di \emph{post-processing} permettono all'utente di salvare e visualizzare la soluzione del problema parabolico.
Nello specifico, il metodo \texttt{printSol} effettua semplicemente una chiamata al metodo \texttt{save\_to\_file} della classe \texttt{Mesh}:
\begin{c++}
template<class BIL,class F>
void Problem<BIL,F>::printSol(char const * filename) {
    mesh.save_to_file(_u, filename);
}
\end{c++}
In questo modo viene creato un file di testo che memorizza le coordinate dei nodi, oltre ai corrispondenti valori della soluzione per ogni passo temporale.\\
A partire da questo file, in un secondo momento, è poi possibile ottenere le immagini della soluzione nei singoli istanti temporali, oppure un filmato che rappresenti l'evoluzione nel tempo della soluzione stessa, grazie al supporto del già citato programma di visualizzazione \texttt{gnuplot}.
I metodi che si occupano di generare gli script necessari, atti a svolgere queste operazioni, sono rispettivamente \texttt{printGnuplotScript} e \texttt{printAviScript}.
\begin{itemize}
	\item Il metodo \texttt{printGnuplotScript} genera uno script per la bash che si occupa di lanciare il software \texttt{gnuplot} e di eseguire i comandi necessari alla visualizzazione grafica della soluzione in un certo passo temporale.
	Tutti i parametri geometrici vengono adattati automaticamente al particolare caso preso in considerazione, grazie al fatto che il metodo, trovandosi all'interno della classe \texttt{Problem}, è in grado di conoscere tutte le informazioni necessarie a generare immagini ottimali.
	\item Il metodo \texttt{printAviScript}, invece, genera uno script per la bash che permette di ottenere un filmato, in formato \texttt{.avi}, dell'evoluzione nel tempo della soluzione.
	In questo caso le singole immagini sono ottenute sempre attraverso l'impiego di \texttt{gnuplot}, mentre la funzione di \emph{encoding} è svolta dal software gratuito \texttt{MEncoder}, il quale comprende anche lo strumento di riproduzione \texttt{MPlayer}.
\end{itemize}
Grazie a queste funzionalità, l'utente ha dunque a disposizione degli strumenti che si adattano alle caratteristiche della soluzione ottenuta e che, allo stesso tempo, forniscono la possibilità di essere modificati secondo le necessità del problema specifico che si intende risolvere.
Ricordiamo, infine, che per rendere eseguibili i diversi script è prima necessario modificare i permessi dei file, per esempio attraverso il comando della bash \texttt{chmod}.


\section{Strumenti software e librerie esterne}

Come accennavamo all'inizio del capitolo, durante la fase di sviluppo della libreria \texttt{Heat-Spectral} ci siamo avvalsi dell'aiuto di altri supporti software.
In particolare abbiamo introdotto il sistema di controllo di versione \texttt{subversion}, il generatore automatico di documentazione \texttt{doxygen}, la libreria per maneggiare le matrici sparse \texttt{UMFPACK}, oltre agli strumenti di visualizzazione ed encoding \texttt{gnuplot} e \texttt{MEncoder}.
A questo punto passiamo ad analizzare in dettaglio ciascuno di questi supporti, mostrando in che modo hanno contribuito alla realizzazione della nostra libreria.

\subsection{Subversion}

Il sistema di controllo di versione \texttt{subversion} è un programma multi-piattaforma distribuito con licenza di software libero \emph{Apache}.
Il meccanismo di funzionamento di questa applicazione è relativamente semplice.
Di fatto, la versione definitiva del progetto è conservata in un \emph{repository} centrale, dal quale tutti gli sviluppatori possono ottenerne una copia di lavoro, sulla quale apporre eventuali modifiche.
In un secondo momento, una volta completata una parte del lavoro, è possibile caricare le modifiche sul server, dove sarà il sistema di controllo a occuparsi di sincronizzare i contributi di ogni sviluppatore.
Le caratteristiche di un sistema costruito in questo modo sono essenzialmente tre:
\begin{itemize}
	\item il vantaggio, da parte dei vari programmatori, di gestire in maniera simultanea il software;
	\item la possibilità di mantenere traccia di ogni modifica attraverso la memorizzazione di tutte le versioni;
	\item la sincronizzazione e l'aggiornamento costante del lavoro in corso d'opera.
\end{itemize}
Per quanto riguarda il primo punto, è indispensabile che il sistema consenta a più sviluppatori di lavorare contemporaneamente allo stesso progetto.
In questo modo, tutti gli eventuali problemi che si possono generare quando più persone lavorano allo stesso codice vengono trattati in maniera automatica.
In pratica è lo stesso \texttt{subversion} che identifica i conflitti, oltre a fornire all'utente gli strumenti necessari a risolverli.
La caratteristica fondamentale è che nulla viene perduto o sovrascritto, ma sono gli sviluppatori a decidere quale versione mantenere tra quelle concorrenti.\\
D'altra parte, il fatto di avere un archivio di tutte le modifiche è un grosso vantaggio; infatti consente sempre di tornare a una vecchia versione stabile, nel caso in cui si siano apportati degli errori che hanno danneggiato il codice.
Inoltre, il sistema tiene traccia dell'autore che compie le modifiche, oltre a consentire un'immediata comparazione tra diverse versioni dello stesso file.
In questo modo si possono tenere sotto controllo i progressi svolti e si possono identificare rapidamente ogni modifica apportata.\\
Infine, il fatto di poter sincronizzare costantemente il lavoro consente ai vari sviluppatori di rimanere aggiornati sulle modifiche compiute da altri, oltre a rendere tutti quanti partecipi del proprio lavoro.\\
Non abbiamo intenzione di scendere nei dettagli relativi alla sintassi del client ma, nel caso si voglia approfondire l'argomento, è possibile consultare il manuale ufficiale in~\cite{SUBVERSION}.

\subsection{Doxygen}

Il generatore automatico di documentazione \texttt{doxygen} è un progetto \emph{open source}, distribuito sotto licenza GPL, che opera con tutti i principali linguaggi di programmazione.
Come dicevamo nell'introduzione, è uno strumento che si rivela molto utile per risolvere l'annoso problema del manuale d'uso per il codice a cui si sta lavorando.
Grazie alle sue caratteristiche, infatti, è in grado di produrre una ricca documentazione, semplicemente partendo dai commenti inseriti all'interno dei sorgenti e, naturalmente, sfruttando le strutture dati che compongono il nostro progetto.
La documentazione può essere generata in formato HTML, per la navigazione mediante un \emph{browser} qualsiasi, oppure in formato \LaTeX, in modo da ottenere un manuale stampabile.
In effetti, la complessità di \texttt{doxygen} gli consente di supportare anche numerosi altri formati, tra cui PDF, PostScript, RTF e \emph{man page} di Unix.\\
Dicevamo che le descrizioni relative alle singole classi e a ciascun metodo sono inserite direttamente nei file sorgenti come commenti.
È però necessario rispettare alcuni accorgimenti di sintassi richiesti da \texttt{doxygen}, per fare in modo che la procedura di produzione del manuale vada a buon fine.
Inoltre, l'applicazione è in grado di generare in maniera automatica anche i diagrammi delle dipendenze tra le varie strutture, in modo da rendere meno onerosa la comprensione del funzionamento del programma che vuole documentare.
A ogni progetto, \texttt{doxygen} associa un file di configurazione, che contiene tutte le impostazioni necessarie a generare il manuale.
In questo modo, attribuendo particolari valori ai parametri richiesti, come il nome del progetto oppure la posizione dei sorgenti, la documentazione può essere prodotta in maniera automatica.
Un file di configurazione di default può essere generato attraverso il comando
\begin{verbatim}
doxygen -g <config-file>
\end{verbatim}
eseguito da riga di comando.
Invece, per quanto riguarda l'esecuzione dell'applicazione, può essere facilmente avviata per mezzo del comando
\begin{verbatim}
doxygen <config-file>
\end{verbatim}
In ogni caso, per tutti i dettagli relativi alla configurazione e alla documentazione di \texttt{doxygen}, si rimanda al sito ufficiale del progetto\footnote{\url{http://www.stack.nl/~dimitri/doxygen/}}.

\subsection{UMFPACK}

La libreria \texttt{UMFPACK} è un insieme di metodi che permette di risolvere sistemi lineari del tipo
\[
	\A\vec{x}=\vec{b}
\]
nel caso in cui $\A$ sia una matrice sparsa e non simmetrica.
L'insieme di strumenti contenuti in questo pacchetto è in grado di effettuare una fattorizzazione LU della matrice, che in realtà può essere sia quadrata sia rettangolare, singolare oppure non singolare, composta di numeri reali oltre che complessi.\\
La libreria è composta da $32$ routine che possono essere richiamate direttamente dall'utente, oltre che dal file \texttt{umfpack.h}, che deve essere incluso nel programma in linguaggio C++ che utilizza la libreria.
Nello specifico, per ciascuna routine sono disponibili $4$ versioni, ognuna per un tipo diverso: valori interi o \emph{floating-point}, reali oppure complessi.\\
Oltre a questo, esistono $5$ routine principali che intendiamo descrivere in questa breve introduzione:
\begin{itemize}
	\item \texttt{umfpack\_*\_symbolic}. Effettua un pre-riordinamento delle colonne di $\A$, in modo da ridurne il riempimento.
	Dopodiché ritorna un oggetto di tipo \texttt{Symbolic} come un puntatore a \texttt{void}.
	Questo oggetto contiene l'analisi simbolica ed è necessario per la fattorizzazione numerica.
	Di fatto determina anche i valori limite degli elementi non nulli di $\mathrm{L}$ e $\mathrm{U}$, oltre alle operazioni necessarie e la richiesta di memoria per il successivo metodo \texttt{umfpack\_*\_numeric}.
	Questa routine richiede uno spazio in memoria naturalmente proporzionale al numero di elementi non nulli di $\A$.
	\item \texttt{umfpack\_*\_numeric}. Fattorizza la matrice del tipo $\mathrm{PRAQ}$ nel prodotto $\mathrm{LU}$, dove $\mathrm{P}$ e $\mathrm{Q}$ sono matrici di permutazione, $\mathrm{R}$ è una matrice diagonale di fattori di scala, $\mathrm{L}$ è una matrice triangolare inferiore con diagonale unitaria e $\mathrm{U}$ è una matrice triangolare superiore.
	Il metodo sfrutta l'analisi simbolica effettuata dalla routine \texttt{umfpack\_*\_symbolic} per restituire un oggetto di tipo \texttt{Numeric} come puntatore a \texttt{void}.
	Questo oggetto contiene la fattorizzazione LU e viene usato dalla routine \texttt{umfpack\_*\_solve}.
	\item \texttt{umfpack\_*\_solve}. Risolve il sistema lineare $\A\vec{x}=\vec{b}$ utilizzando la fattorizzazione numerica effettuata da \texttt{umfpack\_*\_numeric}, dove la matrice $\A$ deve essere quadrata e non singolare.
	\item \texttt{umfpack\_*\_free\_symbolic}. Semplicemente libera lo spazio di memoria occupato dall'oggetto \texttt{Symbolic}.
	\item \texttt{umfpack\_*\_free\_numeric}. Come nel caso precedente, viene utilizzata per liberare spazio di memoria, ma in questo caso occupato dall'oggetto \texttt{Numeric}.
\end{itemize}
Come ampiamente discusso in precedenza, tutte le matrici che entrano in gioco nella risoluzione algebrica del problema differenziale sono memorizzate nel formato sparso \texttt{CSC}, e dunque utilizzano le strutture disponibili in questa libreria.
In ogni caso, per tutti i dettagli relativi a questo tipo di formato, rimandiamo a \cite{FORMAGGIA}.

\subsection{Gnuplot}

Il programma a linea di comando \texttt{gnuplot} permette di rappresentare graficamente funzioni espresse per via analitica, ma anche file di dati ottenuti come risultati di codici numerici.
Il progetto è naturalmente multi-piattaforma, ed è distribuito con una licenza \emph{Creative Commons} che ne determina la natura di software libero.
Sono infatti queste due caratteristiche, assolutamente irrinunciabili quando si lavora alla programmazione di tipo scientifico, che hanno determinato la scelta di questo supporto piuttosto che di un altro.
Inoltre, \texttt{gnuplot} è in grado di realizzare grafici sia bidimensionali sia tridimensionali, oltre ad avere la possibilità di rappresentare punti, linee, superfici, campi vettoriali e stringhe di testo.
Per quanto riguarda il supporto, sul quale tutte queste rappresentazioni possono essere visualizzate, esistono diverse possibilità.
È permesso a \texttt{gnuplot} di interagire direttamente con lo schermo della macchina sulla quale si sta lavorando, oppure può stampare le immagini direttamente su terminali esterni o remoti.
Infine, vi è anche la possibilità di salvare i grafici su file, sfruttando i numerosi formati a disposizione, dai più classici (jpeg o png) a quelli vettoriali (eps, pdf o svg).\\
Come dicevamo, la scelta dello strumento di interfaccia grafica è ricaduta su \texttt{gnuplot} proprio a causa della sua estrema flessibilità.
Essendo un programma a linea di comando, è stato possibile inserire le istruzioni direttamente nei sorgenti della nostra libreria, in modo da generare uno script di facile utilizzo per l'utente finale.
Tale script è chiaramente modificabile in un secondo tempo, così da sfruttare ogni caratteristica del software, per adattarsi alle esigenze dell'utente.\\
Passiamo ora a illustrare in che modo è stato implementato lo script che permette a \texttt{gnuplot} di generare i grafici delle soluzioni.
Innanzitutto è necessario scegliere un tipo di terminale per l'output.
Nel nostro caso abbiamo individuato nel terminale \texttt{pdf} la soluzione ottimale, poiché permette di ottenere immagini in formato vettoriale, utilizzabili direttamente da \texttt{pdflatex} per la compilazione di relazioni o presentazioni.
Tutte le istruzioni che precedono il comando di \emph{plot}, e dunque che consentono di impostare tutte le caratteristiche che l'immagine dovrà avere alla fine, vengono definite attraverso il comando \texttt{set}.
Ad esempio, una tipica configurazione potrebbe essere la seguente:
\begin{verbatim}
set terminal pdf enhanced color fsize 8 size 5,3
set output 'filename.pdf'
set title 'title'
set xrange [0:6]
set yrange [0:6]
set zrange [0:1]
set palette rgbformulae 21,22,23
\end{verbatim}
Queste istruzioni si occupano, rispettivamente, di definire il tipo di terminale con tutte le sue caratteristiche, il nome del file di output, il titolo dell'immagine, gli intervalli in cui rappresentare le tre variabili spaziali e la scala dei colori necessari a rappresentare l'entità della soluzione.
In realtà, per poter accoppiare questo tipo di interfaccia con la nostra libreria, è stato necessario sfruttare le variabili interne del nostro codice, in modo da generare intervalli corretti di visualizzazione, oltre a una scala di colore coerente con la soluzione numerica.
A questo punto, a seconda dell'obiettivo che si vuole ottenere, esiste la possibilità di scegliere se plottare un'immagine bidimensionale oppure tridimensionale, rispettivamente mediante il comando \texttt{set pm3d map} o \texttt{set pm3d at s}.
Infine si procede alla vera e propria generazione dell'immagine, attraverso il comando \texttt{splot}, che si occupa di rappresentare dati in più dimensioni.\\
Questa descrizione sommaria vuole essere soltanto un'indicazione per l'utilizzo delle funzionalità di \texttt{gnuplot}.
La versione completa del metodo che si occupa di generare lo script è stata ampiamente documentata nel manuale della libreria, oltre ad aver accennato a una sua descrizione nelle precedenti sezioni.

\subsection{MEncoder}

Secondo quando illustrato nelle precedenti sezioni, oltre a un'interfaccia grafica, abbiamo dotato la nostra libreria della possibilità di generare in automatico dei filmati, rappresentanti l'evolversi della soluzione nel tempo.
A tal proposito abbiamo pensato che lo strumento ideale, che potesse fare al nostro caso, fosse \texttt{MEncoder}, un \emph{encoder} video multi-piattaforma che permette di generare filmati nei formati più disparati.
Incidentalmente, il pacchetto del software contiene anche un riproduttore video, \texttt{MPlayer}, in grado di sfruttare la tecnologia dei \emph{codec} più moderni.
Entrambi i programmi hanno l'enorme vantaggio di essere distribuiti sotto licenza GNU, ed è anche per questo che sono stati scelti per interfacciarsi con la nostra libreria.
Inoltre, hanno la caratteristica di poter essere eseguiti su diverse architetture e sulle piattaforme più disparate, ma sempre con la stessa efficacia e funzionalità.
Sono queste le ragioni che ne determinano il primato su tutti gli altri \emph{encoder} nel panorama delle applicazioni \emph{open source}.\\
Analogamente al caso di \texttt{gnuplot}, le istruzioni di \emph{encoding} non vengono eseguite direttamente dal nostro codice.
In realtà viene generato uno script, da una funzione della nostra libreria, che poi può essere eseguito in un secondo momento dall'utente.
Inoltre, con \texttt{gnuplot} è possibile ottenere un'immagine della soluzione, riferita solamente a un ben preciso istante temporale.
Invece, attraverso il filmato che \texttt{MEncoder} può costruire, è possibile visualizzare in sequenza tutte queste immagini, in modo da ottenere una rappresentazione dell'evolversi del problema parabolico di forte impatto visivo.\\
A questo punto, l'unica operazione da svolgere è salvare in formato \texttt{jpg} le immagini relative a ogni singolo istante temporale, generate grazie a \texttt{gnuplot}.
Poi, con una semplice istruzione si ottiene il filmato:
\begin{verbatim}
mencoder 'mf://*.jpg' -mf fps=25 -o movie.avi 
         -ovc lavc -lavcopts vcodec=mpeg4:vbitrate=800
\end{verbatim}
Il fatto di impostare alcuni parametri permette di scegliere tutte le caratteristiche principali che il filmato in fase di creazione dovrà avere.
In particolare, relativamente al nostro obiettivo, abbiamo selezionato il numero di \emph{frame} per secondo, il tipo di \emph{codec} e il \emph{bitrate} del video.
In ogni caso, tutte le informazioni necessarie a comprendere il modo in cui l'\emph{encoder} si interfaccia con la nostra libreria si possono trovare nella documentazione.
Ricordiamo infatti che lo script, capace di eseguire l'operazione di \emph{encoding}, è generato in maniera automatica dal metodo \texttt{printAviScript}.
